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We propose a new simulation computational method to solve the reduced BCS Hamiltonian based 
on spin analogy and submatrix diagonalization. Then we further apply this method to solve su- 
perconducting energy gap and the results are well consistent with those obtained by Bogoliubov 
transformation method. The exponential problem of 2'^-dimension matrix is reduced to the poly- 
nomial problem of A''-dimension matrix. It is essential to validate this method on a real quantum 
computer and is helpful to understand the many-body quantum theory. 



PACS numbers: 74.20.Fg, 03.67.Lx 

■ BCS theorjii and its subsequent extension is a well established theory to explain the mechanism of supercon- 
I ducting property. With two gross simplifications: the free electron approximation and the effective interaction 

approximation^.^, a simplified BCS model is obtained and described by the reduced BCS Hamiltonian. There has 
been much work on solving this Hamiltonian. The mean field method is exact in the limit of large number of electrons 
where fluctuation can be neglected but disabled in the case of small number of electrons. Since Richardson's work4 
in the 60's to now, the exactly solvable BCS Hamiltonian attracts much attention in connection with the problems in 
lO different areas of physics such as superconductivity, nuclear physics, physics of ultrasmall metallic grains. 

■ Recently in L.-A. Wu et al's paper— an NMR experiment scheme performing a polynomial-time simulation of pair- 
ing model was reported. Based on this work we propose an explicit theory method to diagonalize the reduced BCS 
Hamiltonian through the spin analogy and submatrix diagonalization. Compared with the conventional method it is 
more useful in solving practical problem. The problem is solved in the spin space, which is convenient related to the 

, qubit system. It gives a senseful alive method, quantum simulation, instead of the numerical diagonalization calcula- 
tion. And it shows the potential to solve many-body problem by quantum computer. In fact more and more people 
C ' concentrate on the research of simulating other physics systems by quantum computer^. The experimental quantum 
• simulations about quantum harmonic oscillator^, three-spin artifical Hamiltonian? and migration of excitation in a 
' one-dimensional chain^ et al. have been realized. Recently a relative experiment is performed to get the eigenvalues 
of the BCS Hamiltonian through selecting a proper initial state and realizing Hamiltonian evolutioniiii'^ 
The exact solvable model, i.e. the reduced BCS Hamiltonian considered in this paper is^i^iiS: 
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where n±m = Cj_„c±m are the electron number operators, c}„(cm) is the fermionic creation (annihilation) operator. 
The coupling coefficient is simplified as a constant V~-^. Note that the summation indexes to = 1, 2, ■ • • ,N represent 
^ . , all of relevant quantum numbers, and the electron pairs are labelled by the the quantum number m and —to, according 
' to the Cooper pair situation where the paired electrons have equal energies but opposite momenta and spins: to = 
( fc , t ) and —TO = (— fc , J,). Introduce the pair creation operator = cl^cL^n and the pair annihilation operator 
b,n = c-mCm- So One can write the Hamiltonian as^^: 

N ^ N 

HbCS - E T^"" + ""-"'^ ^^rnh (2) 

m=l m,l=l 

where = Em ~ £f is the free electron kinetic energy from Fermi surface (e f is the Fermi energy) . There are two 
possible cases for every pair state m: "occupation" and "empty" , which are denoted respectively by: 

where the spin up state xi indicates "occupation" and the spin down state xo indicates "empty" . Obviously, — {a^ — 
icry)xi — Xo and 2^<^^ ~ ^'^y)Xo = 0, then we can get the so-called spin-analogy corresponding of the pair annihilation 
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operator bm as 
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In the same way, the spin-analogy corresponding of the pair creation operator bl^ becomes 

=C7+ (5) 

From the pair number operator nm + n^rn has the eigenvalue 2 (which represents the electron number in every Cooper 
pair) when operating on xi, and when operating on xoi it follows that 

+ n.m E!) = 1 + ^i"^ (6) 

In fact, the fermionic pair operators satisfy the commutation algebra: sl{2) — {fomj&Li^m + ?^-)ri — l}, J'-C- sl[2) = 
{cr~ , (7^, cr^}. From formulas Q)-© one can express H^cs in terms of the spin operators: 
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where em = -V. 

In fact the spin analogy of the BCS Hamiltonian is well known and exact diagonalization of the pairing model in 
the spin space has been carried out in several previous works^^. In this paper we propose a computational simulation 
method which is potential to realize in future with the development of quantum computer. The primary advantage 
of our method lies in the practical realization in experiments. Especially we can solve superconducting energy gap 
by this simulation method conveniently as following paragraphs. It is more practicably than other solution of energy 
gap, because it can simplify a 2^-dimcnsion problem to an iV-dimension problem. We know eigenvalues may not 
be solvable for high dimension matrix in principle. Now a 2^-dimension problem, exponential problem(EP) can be 
simplified an iV-dimension problem, polynomial problcm(PP). In the following part we will describe how to transform 
EP to PP in detail. 

Firstly the total Hamiltonian iJspin will be expressed as the direct-sum of a set of submatrices^. 

^^spin = HsubO © Hsubl © Hsub2 © ' ' ' ® HgubN (7) 

The system states with the same spin- up state number form an absolute subspace. The subscripts subO, subl, ■ ■ ■ subN 
representing the number of the spin- up state in the corresponding subspace are respectively 0, 1, • ■ • , A^. Secondly we 
will prove that the eigenvalues of Hsubi in subl submatrix justly are the eigenvalues of -ffecs- 

Note the 2^ x 2^ operator as 

h{N, m) = ® (^J 0^ ® i<»{N-m) (to = 1^ 2, • • ■ , iV) (8) 

whose i-th diagonal element is noted as h{N^m)[i] {i — 1,2, .....N). It is easy to see the non-diagonal elements of 
h{N, m) are zero. 

Lemma 1 For the Hamiltonian as h{N,ra) = I®!™-!) ® (^q 8) I^C^-™), the (2^ - 2^-*)-t/i element satisfies 

Wm)[2^-2^-''] =(5,™ (i,m= 1,2,.-. ,iV) (9) 
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Proof. We will prove this lemma by mathematic induction. 

When = 3, it is easy to get three matrices, m value is 1, 2, 3 respectively. 
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Here the non-diagonal elements are all 0. These diagonal elements with ' ' are h{3, m)[2^ — 2^ It is easy to validate 
h{3, m)[2^ — 2^~'] = and the last diagonal element /i(3, m)[2'^] is zero. If when N = L, 



is right and 



h{L, m) [2^ - 2^-'] =S^m (z, m = 1, 2, • • • , L) 



h{L,m)[2^] = 



(10) 



(11) 



then we should examine whether h{L + 1, m)[2^^-^ — 2^+^ — Sim {i,m — 1,2, L + 1) is right when N ^ L + 1. 
We will discuss it in two cases: m ^ L and m ^ L + 1. 
(1) m s; L : 



h{L + l,m) = «> (^J ® 

= h{L,m) (g) I 
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From Ea. (|12|l it is easy to see 



/i(L + 1,to)[2x] = h{L + l,m)[2x-l] ^h{L + l,m)[ 



(12) 
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/i(L + 1, m)[2^+^ - 2^+^-'] = h{L, m)[2^ - 2^"^] = (5„„ (i = 1, 2, • • • , L) 

Next we should also know the value of /i(L + l,m)[2^+^-l]. From Ea. (fTT|l andEq.(0 there is /i(L-|-l, m)[2^+^ - 1] = 
h{L,m)[2^] = 0. From above discussion for m < L{N = L + 1), the equality h{L + l,m)[2^+i - 2^+^"*] = 4„i is 
valid. 

(2) m = L + 1: 



/i(L + l,L + l) = 1®^® 
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From the above expressions the (2^+^ — l)-th diagonal element h{L + 1, L + 1)[2^+^ — 1] = 1, that is to say 

h{L + 1,L + l)[2^+i - 2^+^-'] = l{i = L + l) 

And then for z L + 1, 

h{L + 1,L + l)[2^+i - 2^+1-*] = 
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because the even-th diagonal elements are all zero obviously. So for m = L+1 there is also h{L+l, m)[2^ — 2-^"*] = 5i,n- 
From the discussion (1) and (2) we have proved that when N = L + 1, h{L + l,m)[2^+^ — 2^+^^'] = (5^^ (i,m = 
1, 2, • ■ • , L + 1) is valid. So for ViV ^ 3 [N is the natural number), there is the equality 

h{N, m) [2^ - 2^-'] =5,^ (z, m = 1, 2, • ■ • , TV) 



After the preparation we will prove that the eigenvalues of Hgubi in subl subspace justly are the energy spectrum of 
quasiparticle excitation of ffecs- Firstly set a diagonal Hamiltonian as iJ^'^'s = i J2m=i EmijmJm + J-ml-m) , 7m 7™ 
and 7l,„7-m are the quasiparticle number operators. According to the previous analogy rule of number operators 
Eq.® the spin-analogy form of H'^^^^ is 

N 
m — 1 

-^spTn's submatrix in subl subspace is denoted as -ff^^^f and the i-th diagonal element of -ff^^^^i ^ ^suwl*]- 
can find the i-th diagonal element of H^^f is the (2-^ — 2^~*)-th diagonal element of the total Hamiltonian H^p^^- 
According to Eq.® 

fjdiagr-l _ rrdiagrrjTV nN—i] 
-"siifell'J — -"spin ^ J 

N 

m=l 
N 

= E, 

So the spin-analogy Hamiltonian of a diagonal BCS Hamiltonian i/'^'^s subl subspace has the same eigenvalues 
as the diagonal BCS Hamiltonian's energy spectrum. We can deduce further that the eigenvalues of Hsubi are justly 
the eigenvalues of Hbcs whether which is diagonal or not, because Hbcs can be written as the diagonal form like 

ffdiag generally, i.e. Hbcs — > 5 Sm=i ^milmlm + l-m'y-m) H , while we need not care about how to obtain the 

diagonal form Hamiltonian. Consequently it implies that if we have diagonalizcd Hsubi we can get the energy spectrum 
of quasiparticle excitation of Hbcs and get energy gap further avoiding complex computation. One of the classical 
methods to diagonalize Hbcs, Bogoliubov transformation method is available under the mean-field approximation, 
which is not exact especially in the case of limited N. Another famous method about the exact solution of BCS 
Hamiltonian has been proposed in the 60 's by Richardson^. He considered the system with M ^ N pairs electrons 
and constructed a set of operators Rj (j ~ 1, • • ■ , -/V) commuted with Hbcs, finally gave the expression of the 
eigenvalues Aj of Rj through solving the M coupled algebraic equations. Aj is not yet the eigenvalue of Hbcs- We 
can more directly and simply give the energy spectrum of Hbcs- 

Thus it can be seen the idea of diagonalizing Hsubi instead of solving the eigenvalues of Hbcs directly is better than 
those classical ones. Now a key problem that how to get the submatrix Hs^bi is placed to the front. The correlative 
workiSi in our group has proved that the general form of Hsubi ■ 

Hsubi[i,i]^e,, Hsubi[iJ] = -V = 1,2,--- ,N; ij^j) (14) 

Hsubiihj] is the matrix element of Hsubi- 

Finally it is necessary to check the methods in numerical computation. Here we will compare our solution with 
the result of the mean-field approximation by the value of superconducting energy gap A (T = OK). According to 
the physics meaning of A, the energy required to excite at least a quasiparticle from the Fermi surface, the energy of 
the element excitation is written as = {^m + A-^)^/^i^ In fact the element excitation energy is also the eigenvalue 
of Hbcs- After getting the eigenvalue of Hsubi, i-S- the eigenvalue of Hbcs, we can get the value of A by solving 
equation Em = (Cm + A^)^/^. But in order to get rid of the effect of energy zero the equation 



(15) 
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is used to solve A in practice, because the difference of eigenvalues dosen't depend on the energy zero and that we 
find the energy difference between the ground and the first excited state E2 — Ei is by far larger than — E2, E4 — 
i?3, • • • , En — En-1 in the course of the numerical computation. Here ^1 < S.2 < ■ ■ ■ < S,n and Ei < E2 < ■ ■ ■ < En. 



The another kind of solution used in comparing is the following energy gap equationi^: 

1 



(16) 



Here we consider the reduced BCS model whose energies are given for simplicity by = mdi^^^, here S is 
the average level spacing which is inversely proportional to the size of the grains. In the strong coupling regime, 
corresponding to large grains or strong coupling constants, (5 <C A. In the weak coupling region, corresponding to 
small grains or small coupling constants, S 3> A^^. From much research about ultrasmall superconducting grains^S, 
the mean-field theory is not suitable in the weak coupling region. It has been proved that the corrections to the 
mean- field results are small in large grains become important in the opposite limit So in this paper we carry out 
the numerical computation in the first case, (5 <C A. We also take the coupling constants V — X6 in order to discuss 
conveniently. In order to give the numerical pictures, we suppose the value of V by the rough estimate. From BCS 
theory, the Cooper pair lies in the attraction area, i.e. < < Hzuef^- For metal Debye energy Kwd ~ 10~^eV, we 
can set V ^ 2x 10~^eV by rough estimate. The estimate process is put to the later appendix. Another two variables 
A and N are taken as the independent variables of the energy gap. 

We list our results in diagrams. In FIG^ setting A = 10, the energy gap is plotted as the function of energy level 
number, which is the mono-increasing function of the energy level number N . 
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FIG. 1: The energy gap A as a function of the number of energy level, A'^ . We give the comparison between two results by the 
different methods. The solid line is the solution of the energy-gap equation; the dashed line is our result by spin analogy and 
diagonalizing submatrix. Here we choose V ~ 2 X lO^^eV, A = 10 . 

The relative error between our result and that of the Ea. (|16|l is not more than 5% in the range from = 2 to iV = 100 
and fixed A(= 10). It shows that the result from our method is well consistent with the solution of the energy gap 
equation. In order to check the universality of this new method, we also give the dependence relation between the 
energy gap and the level spacing, see FIG [3 

It is clear to see A changes gently with A when A is large enough. That is to say, A is almost independent of 5 when 
5 is small enough. It shows the rationality of 5 ^ A on the inverse hand. Obviously, in FIG|2 when A is larger than 
80, the relative error is less than 1.1%. 

We also consider a small departure from the fermi surface, that is = ^0 -f m5. is a small value, ^0 5. In the 
following discussion we note ^0 = b is the natural number. The small departure from the fermi surface reduces the 
energy gap and energy gap is not a real root when the departure reaches a critical value, see FIG|31 

According to above comparison we know two results are consistent well, while our method to solve the energy 
spectrum doesn't include approximation, which indicates that our result includes that obtained by mean field theory 
and is superior to it. 

In summary, we have proposed an exact numerical simulation method to calculate the energy spectrum of the 
reduced BCS Hamiltonian by spin analogy and diagonalizing submatrix. A numerical computation to verify the 
validity of our computational method is given. We make a comparison between our method and energy gap equation 
Ea. (|16|l . and two results are well consistent in numerical computation. By examining the change of the energy gap 
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FIG. 2: When A'^ = 20, Energy gap A is plotted as the function of A. We give the comparison between two resuhs: the sohd 
hne is the solution of the energy-gap equation; the dashed line is our result by spin analogy and diagonalizing submatrix. Here 
we choose F ~ 2 x 10"'' eV. 
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FIG. 3: When 1/ ~ 2 x 10~ eV, A — 10, A'^ = 10 the energy gap A is ploted as the function of b through the quantum 
simulation. When b > 56, there is not the real root for A. The result is obtained by the quantum simulation. Similar result 
will be obtained by the energy gap equation Eg. 1161 . 

value under the change of the parameter, we include the excellent consistency between the two results by the different 
methods is independent on the particular parameter. It implies that one can implement this quantum simulation on 
a quantum computer and the result will be believable. Currently a new experiment about 2-qubit simulation of the 
pairing Hamiltonian on an NMR quantum computer has been realized and get the energy spectrum of the pairing 
Hamiltonian successfuUjiifi. With the development of quantum computer, especially the manipulation and control of 
multi-qubit system, this new simulation computation method has the great potential in practical application. 
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APPENDIX A: APPENDIX: ESTIMATE ABOUT V 



For metal element g{0)V « 0.2 ^ 0.2^^, 17(0) is the state density which have some spin directions on the Fermi 
surface. 



/"OO 

/ 2g{e)d\ = 1 
^0 



(Al) 
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According to the assumption 



g{e) = 0, if e > hwo 
g{e) — g{0), if e < hruo 



formalism (^3J can be written as 



/ 2g{Q)dh = 1 

^0 



so from BCS theory, the Cooper pair hes in the attraction area, ie < ^fc < hwu- On substitution of Hwd ~ 10 ^eV— 
we can estimate g(0) « 10^, so y ~ 2 x lO^^eV. 
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